﻿/* Copyright 2008 dnAnalytics Project.
 *
 * Contributors to this file:
 * Jurgen Van Gael
 * 
 * The implementation of the incomplete gamma function is based on Numerical Recipes 3d edition.
 *
 * This file is part of dnAnalytics.  dnAnalytics is licensed under the 
 * Microsoft Public License. See License.txt for a complete copy of the
 * license.
 */

using System;

namespace dnAnalytics.Math
{
    public static partial class SpecialFunctions
    {
        // Decides when to evaluate using quadrature methods.
        private const double mIncGammaSwitch = 100.0;

        // Quadrature parameters.

        // Quadrature parameters.
        private static readonly double[] mIncGammaW =
            {
                0.0055657196642445571, 0.012915947284065419,
                0.020181515297735382, 0.027298621498568734,
                0.034213810770299537, 0.040875750923643261,
                0.047235083490265582, 0.053244713977759692,
                0.058860144245324798, 0.064039797355015485,
                0.068745323835736408, 0.072941885005653087,
                0.076598410645870640, 0.079687828912071670,
                0.082187266704339706, 0.084078218979661945,
                0.085346685739338721, 0.085983275670394821
            };

        private static readonly double[] mIncGammaY =
            {
                0.0021695375159141994, 0.011413521097787704,
                0.027972308950302116, 0.051727015600492421,
                0.082502225484340941, 0.12007019910960293,
                0.16415283300752470, 0.21442376986779355,
                0.27051082840644336, 0.33199876341447887,
                0.39843234186401943, 0.46931971407375483,
                0.54413605556657973, 0.62232745288031077,
                0.70331500465597174, 0.78649910768313447,
                0.87126389619061517, 0.95698180152629142
            };

        /// <summary>
        /// Computes the incomplete gamma function P(a,x) using its series representation.
        /// </summary>
        private static double IncompleteGammaSeries(double a, double x, Boolean normalize)
        {
            double sum = 1.0/a;
            double del = 1.0 / a;

            for (double ap = a + 1.0;; ap += 1.0)
            {
                // Save the last sum.
                double prev = sum;

                // Update sum.
                del *= x/ap;
                sum += del;

                // Stop adding terms when adding one doesn't change the representation of sum.
                if (Precision.EqualsWithTolerance(sum, prev, dnAnalytics.Constants.UlpsForComparison))
                {
                    if (normalize)
                    {
                        return sum*System.Math.Exp(-x + a*System.Math.Log(x) - GammaLn(a));
                    }

                    return sum*System.Math.Exp(-x + a*System.Math.Log(x));
                }
            }
        }

        /// <summary>
        /// Computes the incomplete gamma function Q(a,x) using its continued fraction representation.
        /// </summary>
        private static double IncompleteGammaContinuedFraction(double a, double x, Boolean normalize)
        {
            const double FPMIN = 1e-30;

            double b = x + 1.0 - a;
            double c = 1.0/FPMIN;
            double d = 1.0/b;
            double h = d;
            for (int i = 1;; i++)
            {
                double an = -i*(i - a);
                b += 2.0;
                d = an*d + b;

                if (System.Math.Abs(d) < FPMIN)
                {
                    d = FPMIN;
                }
                c = b + an/c;
                if (System.Math.Abs(c) < FPMIN)
                {
                    c = FPMIN;
                }

                d = 1.0/d;
                double del = d*c;
                h *= del;
                if (System.Math.Abs(del - 1.0) <= Double.Epsilon)
                {
                    break;
                }
            }

            if (normalize)
            {
                return System.Math.Exp(-x + a*System.Math.Log(x) - GammaLn(a))*h;
            }

            return System.Math.Exp(-x + a*System.Math.Log(x))*h;
        }

        /// <summary>
        /// Computes the regularized incomplete gamma function using Gaussian Quadrature.
        /// </summary>
        /// <param name="a">The parameter of the gamma function.</param>
        /// <param name="x">The integration upper or lower bound.</param>
        /// <param name="lower">Switch to specify whether to evaluate the lower incomplete gamma or upper incomplete gamma.</param>
        private static double IncompleteGammaQuadrature(double a, double x, Boolean lower)
        {
            double a1 = a - 1.0;
            double lna1 = System.Math.Log(a1);
            double sqrta1 = System.Math.Sqrt(a1);

            double xu;
            if (x > a1)
            {
                xu = System.Math.Max(a1 + 11.5*sqrta1, x + 6.0*sqrta1);
            }
            else
            {
                xu = System.Math.Max(0.0, System.Math.Min(a1 - 7.5*sqrta1, x - 5.0*sqrta1));
            }

            double sum = 0;
            for (int j = 0; j < 18; j++)
            {
                double t = x + (xu - x)*mIncGammaY[j];
                sum += mIncGammaW[j]*System.Math.Exp(-(t - a1) + a1*(System.Math.Log(t) - lna1));
            }

            double ans = sum*(xu - x)*System.Math.Exp(a1*(lna1 - 1.0) - GammaLn(a));
            return (lower ? (ans > 0.0 ? 1.0 - ans : -ans) : (ans >= 0.0 ? ans : 1.0 + ans));
        }

        /// <summary>
        /// Calculates the normalized incomplete gamma function.
        /// </summary>
        /// <param name="a">The parameter of the gamma function.</param>
        /// <param name="x">The integration upper or lower bound.</param>
        /// <param name="lower">Decides whether to evaluate the lower or upper incomplete gamma function.</param>
        public static double IncompleteGamma(double a, double x, Boolean lower)
        {
            if (Double.IsNaN(a) || a <= 0.0)
            {
                throw new ArgumentException("a");
            }
            if (Double.IsNaN(x) || x < 0.0)
            {
                throw new ArgumentException("x");
            }

            if (lower)
            {
                if (x == 0)
                {
                    return 0.0;
                }
                if (Double.IsPositiveInfinity(x))
                {
                    return 1.0;
                }
                if (a >= mIncGammaSwitch)
                {
                    return IncompleteGammaQuadrature(a, x, true);
                }
                if (x < a + 1.0)
                {
                    return IncompleteGammaSeries(a, x, true);
                }

                return 1.0 - IncompleteGammaContinuedFraction(a, x, true);
            }
            if (x == 0)
            {
                return 1.0;
            }
            if (Double.IsPositiveInfinity(x))
            {
                return 0.0;
            }
            if (a >= mIncGammaSwitch)
            {
                return IncompleteGammaQuadrature(a, x, false);
            }
            if (x < a + 1.0)
            {
                return 1.0 - IncompleteGammaSeries(a, x, true);
            }
            return IncompleteGammaContinuedFraction(a, x, true);
        }
    }
}